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Combining renormalization group theoretical ideas with the integral equation approach to fluid 
structure and thermodynamics, the Hierarchical Reference Theory is known to be successful even 
in the vicinity of the critical point and for sub-critical temperatures. We here present a software 
package independent of earlier programs for the application of this theory to simple fluids composed 
of particles interacting via spherically symmetrical pair potentials, restricting ourselves to hard 
sphere reference systems. Using the hard-core Yukawa potential with z — 1.8/a for illustration, we 
discuss our implementation and the results it yields, paying special attention to the core condition 
and emphasizing the decoupling assumption's role. 
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I. INTRODUCTION 



The Hierarchical Reference Theory of fluids (hrt) pioneered by Parola and Reatto has been found well capable 
|^-|llj of describing structural and thermodynamic properties of fluids even in the vicinity of the critical point and 
1 for subcritical temperatures, yielding rigorously flat isotherms in the coexistence region (thus eliminating the need 
for Maxwell constructions) and non-classical values for the critical exponents ||. Still, adoption by a significant 
£^ ■ part of the liquid physics community of this renormalization group (rg) theoretical approach to the integral equation 
description of fluids has largely been lacking so far. While this may be partially attributed to hrt's inherent difficulties 
■ and rather high computational cost, lack of an easy to use yet flexible, well-documented implementation of hrt may 
7— I ' also have played a role. To fill this gap, we have written a program^ suited as a general framework for the exploration 
and application of hrt to simple one-component fluids with hard sphere reference systems with various combinations 
of physical systems, approximations, and solution algorithms. Within the natural limitations of the method, it has 
proved well applicable to a variety of model systems including the hard-core Yukawa (present contribution) and 
hard-core multi-step potentials while most attention has been devoted to the square well fluid jl2| . 

Of course, ours is not the first implementation of hrt for simple one-component fluids: indeed, there has been a 
series of earlier programs j^-Q] by the authors of the theory and their collaborators, but it was the one used in 
henceforth referred to as the "original" implementation, that was a vital step in demonstrating the viability of hrt for 
continuous systems below the critical temperature; though never published or formally released, it has been circulating 
among interested physicists for quite some time. Our new software, on the other hand, differs from its precursors 
in many respects: adoption of a meta language in our version, programming style and documentation-to-code ratio 
may be most obvious, number and nature of hard-coded limitations, important details of the numerical procedure 
and a possible speed gain through generation of customized code might be less apparent. Most importantly, though, 
the original implementation's structure makes experimentation with different combinations of approximations, partial 
differential equation (pde) solving algorithms, parameter settings and physical potentials rather cumbersome; in 
contrast, the fully modular approach adoption of a meta language allowed us to take seems far better suited to a more 
general survey of hrt's numerous attractive features. In addition to the necessary flexibility of our software, great care 
has been taken to ensure the numerical soundness of every step in the calculation and hence of the results produced, 
so that the generation of numerical errors necessarily arising from finite-precision arithmetic is uniformly spread over 
all of the problem's domain. To this end we introduce one central parameter, e#, characterizing the maximum relative 
error introduced at any step; together with a number of criteria relying on e#, this parameter governs virtually all 
of the numerics. Any deviation from this strategy is made explicit, as are all the other approximations entering 
the calculation. Ultimately, our goal was to provide the liquid physics community with a general and versatile yet 
numerically reliable tool for the systematic exploration and assessment of hrt and of the effects introduced by different 
approximations. 



1 Available on the world wide web from ittp: //purl . oclc . org/NET/ar-hrt-1/ 
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This paper is meant to serve a twofold purpose: to present the software we have written and its capabilities, and 
to provide its prospective users with some rudimentary documentation. To this end, after a brief presentation of 
the standard formulation of hrt for one-component fluids and some of the theory's properties as far as they concern 
our implementation (section |d|) , we first give a general outline of our program in section III , only touching upon the 
meta-language it has been written in. Due to our implementation's fully modular design it is only natural to then 
proceed by a discu ssion of t he mo re important of its building blocks and the various approximations they implement 
(sub-sections IIIB through HIE). Presentation and critical assessment of the kind of results that can be attained 
with these and concluding remarks (section IV) are followed by two appendices dealing with some technical aspects 
of the formulation used. 



II. THE THEORY 

While a much more detailed account with additional references can be found in the review article jjj, in what 
follows we want to limit ourselves to only a rough sketch of hrt; in doing so we are going to stress several aspects 
— some of them hardly discussed in the literature — relevant to our implementation of hrt, but only to the extent 
necessary for the discussion thereof. No prior knowledge of hrt is assumed. 

The basic ingredient of hrt, already present in its precursor |Dj ], is the gradual transition from a reference potential 
v (r) to the full potential v(r) = v rcl (r) + w(r) describing the interaction between pairs of particles of a fluid, with 
any one of the intermediate potentials serving as a reference system with respect to which the properties of a successor 
potential are calculated [superscripts always indicate the system a quantity refers to] . In our work we have restricted 
ourselves to the case of a spherically symmetric pure two-body interaction, and we have taken advantage of the 
additional simplifications possible by identifying the reference system with a pure hard sphere system, v lcl — v hs , 
as can always be achieved via the well-known Weeks-Chandler- Andersen scheme G3-C9|. Note however that this 
restriction to a hard sphere reference system is not present in the original implementation of hrt jl 10 ll]]; on the 



other hand, our program's fully modular framework is flexible enough to accommodate any of these extensions should 
the need arise. 

Other than fl3| , hrt achieves the transition from w rof to v in infinitesimally small steps. Inspired by momentum- 
space renormalization group (rg) theory, a cut-off wavenumber Q varying from infinity to zero is introduced, and 
for every Q the potential 

V (Q) _ v rcf + W (Q) jg 

defined such that Fourier components k < Q of the perturbational 
part of the Q-potential are strongly suppressed whereas those for k > Q coincide with those of the original 
potential w. Consequently, the reference system and the fully interacting system are recovered in the limits Q — > 00 
and Q — > 0, respectively: 



v (0) = v _ 



(1) 



The role of the Q-potential just introduced becomes clear when we consider a functional expansion in of ther- 
modynamic and structural properties of the system with pair interaction v^' [a tilde always denotes the Fourier 
transform]: as w^(k < Q) is small, the integrals in the expansion are effectively truncated for k < Q, in keeping 
with the RG picture. 

In principle, the precise manner in which the potential is cut off should not matter, and one can easily conceive of 
many different ways of doing so. On the other hand, for such a procedure to be usable it must not introduce instabilities 
when truncating the hrt hierarchy, which is usually done at the two-particle level. Apart from approaches valid only 
for special types of potentials, we are aware of only two cut-off procedures suitable at least for attractive potentials 
(the standard formulation of hrt can easily be shown to become unstable for w(0) > 0, as implicitly stated already 
in jl]]); in our work we opted for the prescription presented in the review article jl| which seems to have been used 
almost exclusively so far |^J^,^,^|,^7| rather than the smooth cut-off formulation of jjj, the latter being numerically 
cumbersome and predicting non-universal critical exponents. Thus we define the Q-potential 

V (Q) = v rcf + W (Q) by 

*«»(»)-{*«;*>« <2) 

in r-space, w^' differs from w by the addition of a convolution integral, viz. 

(o] M /s 1 f°° ( sin Q Or' - r) sin Q(r' + r)\ , 

w w> (r) = w{r) / r w(r ) dr . 

nr J Q \ r — r r + r J 
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Obviously this Q-potential is a rather artificial function in r-space hardly resembling the full potential except in the 
limits of eq. (|l|); furthermore, the range in r-space over which has to be considered is bound to be much larger 
than that of the original potential — a property immediately carrying over to related quantities, the direct correlation 
functions in particular. As an immediate consequence, numerical Fourier transformations involving the Q-potential 
or any of the correlation functions for the Q-system are computationally expensive and must be treated with extreme 
care; in fact, they should be avoided if possible at all, with obvious repercussions for the implementation of the core 
condition (v. i.). 

In the transition from v lcf to v or, equivalently, from Q = oo to Q = 0, hrt treats the Q-system as reference for the 
properties of the system with infinitesimally lower cut-off Q — dQ; re-summation of terms in the resulting expansions in 
dQ and identification of quantities with a well-defined limit for Q — > finally yields the HRT equations [[y||] : for every 
density g there is a formally exact hierarchy of coupled integro-differential equations involving a suitably modified 
free energy A^ defined as 

egL _ S^L . I ( A0 ) - ^) ( o)) + 1 ( m - *<«(o>) (3 ) 

(cf) = —f3w, f3 = analogously, <j>^ — —flw^Q'), a modified two-particle direct correlation function 

C(Q)( r ) =4 g) (r)+0(r) ~(t> (Q) {r), (4) 

and all higher order correlation functions c$\r), n > 2. The additional terms introduced in A^ and explicitly 
take into account a discontinuity at Q = present in the unmodified free energy Av*> and direct pair correlation 

function cij ; as is apparent from eqs. (||), @) and (|]), modified and unmodified quantities coincide for Q = 0. - 
Furthermore it should be noted that it is customary to include the ideal gas terms in the definition of the Cn '■ for 
the two-particle case this is an additional term of —8(f)/g so that the Ornstein-Zernike Jlgy equation takes the form 

c 2 = -(l/g) - ghc 2 , (5) 

where h(r) is the usual two-particle total correlation function; for the higher order correlation functions cf. Q. 

The full hierarchy the derivation of which we just touched upon yields expressions for the total derivatives with 
respect to Q for C^> and all the c„ , n > 2; for the evolution of A^ we have the particularly simple relation 

dQ\ V J 4tt 2 ^ C(Q)(Q) ) 1 j 

As dc^VdQ, n > 2, involves A iQ) , C (Q) and all higher order correlation functions up to c^ 2 ( so that, in particular, 

dC^/dQ depends on cf ] and cf ] via eq. @)), the equations never decouple and we have to introduce some kind of 
closure. In doing so, it is usually desirable to retain thermodynamic consistency as embodied in the sum rule 



(7) 



rigorously true for the exact solution of the full hierarchy; the derivatives with respect to g present in eq. (Q) then 
mandate the transition from equations at fixed g to a pde in the (Q, p)-plane with boundary conditions supplied 
at two densities, g min and g ma , x . In addition, we need to retain the core condition, viz. g(r) = for r < a where 
g(r) — h(r) + 1 is the pair distribution function; indeed it is one of HRT's main advantages to conserve information on 
all length scales, ranging from the hard-sphere diameter a(g) of the reference system at density g and the associated 
core condition up to the cut-off wave length 1/Q, in the limit Q — > allowing criticality to arise from fluctuations of 
arbitrarily large wave length. 

As noted above, the long-ranged nature of and the correlation functions due to the cutting-off of eq. (||) is a 
strong argument in favor of any closure allowing an approximate implementation of the core-condition without the 
need for costly Fourier transforms. This is a likely reason for the up to now seemingly exclusive use of a closure in 
the spirit of the Lowest- Order ^-ordered Approximation (loga, [ |l9| , p0fl )_or the equivalent Optimized Random-Phase 
Approximation (orpa, ]2l]] ) despite this closure's known deficiencies^ |J : with the argument g silently to be added 
in earlier equations when used within the context of the PDE, we make the ansatz 

CW)(fc, g) = 4>(k, g) + 7 { Q \g) u (k, g) + £M(k, g) , 

W)(k,e) =GW(k,Q)+c? i (k,Q), (8) 
0«>(M =£,?=! 7i Q) (fWn(fc,e), 
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thereby introducing a set of Q-independent basis functions u n and corresponding expansion coefficients 7« . Here, 
u o{ r i q) is chosen proportional to w{r, g) (which has the undesirable effect that, from eqs. (^|) and (|J), the correlation 
function of the Q-system depends upon the full potential w rather than as would be appropriate) and 

normalized so that £io(0, g) — 1; the u n (r, g), n>l, on the other hand, are taken to form a basis for a suitable function 
space over [0, 17(g)]. With these provisions, the problem of implementing both core condition and thermodynamic 

consistency reduces to that of an appropriate choice of the expansion coefficients (g), n > 0, for every point in 
the (Q, p)-plane. With the short-hand notations 

a iQ \g) 

and, for an arbitrary function ip(k, g), 



dQdg 2 \ V ) 



the condition (Q) for thermodynamic consistency is easily re-written as 



dQ VK/ ^ dQ 

and following |(| the core condition can be shown to be equivalent to 

dj ( n Q \g) _ Q 2 4>(Q,g)uj(Q,g) 



(9)-J2^krT l ^(0,g), (10) 



J2muj(k, g)u n (k, g), g] ^-^> = £L — ^ ^ ~ V - ^ L 



n=0 



The latter can be combined with eq. (|T(]) to the more explicit 

E,T=i^ (Q) [uj(k, g) {u n (k, g) - u (k, g) u n (0, g)) , g] a7 fg (g) 

- a^folTW) fu(7c oUnffc nl ol 4- -2l 0(Q,e) ^-(Q.g) • > , I 11 ) 

-a ^jj. LMjlK, «ol«, e;, -I- 27r 2 c«3)(Q, e ) (CCQ>(Q, e )-£(Q, e )) ' J - 

as both the sum rule and the core condition must hold for the reference system, the above evolution equations for the 
7^ are readily supplemented with the initial conditions 7™ f = 0, n > 0. 

As eq. (|ll]) stands, it is no more amenable to direct implementation than the previous formulation; not only must 
this infinite-dimensional matrix equation be truncated to a finite number 1 + N cc of basis functions, but even then 
the resulting integrals need to be evaluated at every Q and g — a tedious process no less demanding than the Fourier 
transformations this approach is meant to replace. What makes this closure manageable, however, is the observation 
that the discontinuity in 1^ 's integrand due to the appearance of instead of the continuous leads to a term 
in di^[ip(k, g), g]/dQ made up of functions evaluated at k — Q alone; following |(|, only this single term is retained, 
leading to the approximation 

d 2 2C(QHQ,g)^(Q,g)~U(Q,g)) 2 

—1^ Mk, g),g] = if>(Q, g) ^ ^- ^ '—3 (12) 

dQ 2n (C(Q)(Q,gj) (c(Q)(Q,e)-mgY) 

leaving out the non-local contribution —2 J2^Lq^-^ [VK^j O) '^n(k 1 g)/c < ^\k 1 g), g](dj^\g) / dQ); with the above ap- 
proximation, the task of evaluating one of the integrals of eq. ( |ll| ) reduces to only an initial integration for the 
reference system followed by the solution of an ordinary differential equation (ode) coupled to the hrt-PDE as well 
as analogous odes for all the other integrals of the Z-type. 

Of course, to fully specify the mathematical problem, the pde must be amended by both initial and boundary 
conditions; while the former take the simple form of vanishing expansion coefficients (v. s.), the latter also impose 
some constraint on the 7n • However, as long as we retain the core condition, such an additional constraint is 

already sufficient to determine the expansion coefficient 7^; unless the 70 so found exactly reproduces eq. JlC|), 
thermodynamic consistency can no longer be imposed without introducing mathematical inconsistencies. By the 
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same token, eq. (11), derived by incorporating the sum rule (0) into the core condition, is no longer valid but must 
be changed to 

EZii {Q) Mk,Q)un(k,g),Q] = -i {Q) [«,-(*, e)fio(*,e),e] 

"^tt 2 C«2)(Q,e)(C<«)(Q, e )-0(Q, e )) ' J - 

to reflect the transition from eq. (Q) to said constraint determining the 7q^ appearing on the above equation's right 
hand side; furthermore, elimination of thermodynamic consistency obviously means decoupling the PDE to a set of 
odes at fixed density. This is exactly what is needed at the boundaries g m in and g max of the density interval considered 
where the p-derivatives defining a^'(g) cannot be evaluated: While specialization to g m in — uniquely determines 
the solution there by virtue of the divergence of the ideal gas term — 1/g in c!> cf (cf. appendix |A|) , some prescription for 
finding 7^, accompanied by the necessary modification of the truncated matrix equation according to eq. (p"3|), must 
be imposed at p max ; as long as we retain eq. jll| ) in the pde's domain's interior and numerically evaluate a^\g) 
there, it is natural to use the loga/orpa condition 7q (g max ) = || at the boundary, which is sufficient to also 

determine dA^\g max )/dQ (or, equivalently, appendix |aJ's f(Q, ftnax)) from the 7« Q ' (ftnax), n>l, alone. 

Unfortunately, it turns out that a scheme retaining c^^(g) in eq. ( fill ) for g m i n < g < g max presents significant 
numerical problems for all but extremely high temperatures, precluding reaching Q = Qq at least for the potentials 
that we have looked at. This is where the so-called "decoupling assumption" comes into play: based upon the different 
ranges of Uo(r) and u n (r), n > 1, the authors of M argue that a^(g) = might be a good approximation, thus 



eliminating the X- integral on the right hand side of eq. (11); it turns out that this change, invariably adopted in 



all later publications, is often sufficient to allow generating a solution all the way to Q = Qo- From our previous 
discussion of conditions imposed on the 7^ it should be obvious that this decoupling assumption is incompatible not 
only with the loga/orpa condition Jq (g) — retained in the original implementation for g = g maK but also with 
thermodynamic consistency (eq. (Q)) altogether; thus we are left with only a few possibilities: we may either retain 
logical consistency by using the decoupling assumption a^\g) — as a closure for the HRT-equations, reducing the 
PDE to a set of odes in Q only; or we may prefer to retain the structure of a PDE so as to make use of thermodynamic 
consistency at least to a certain degree; yet another possibility is to maintain both mathematical and thermodynamic 
consistency by not implementing the core condition at all. The original implementation's approach relying on three 
mutually incompatible concepts, viz. the loga/orpa condition at g max , decoupling, and the compressibility sum rule, 
seems particularly unattractive; at least one should use the decoupling assumption as a boundary condition at high 
density instead. 

Retaining thermodynamic consistency in the form of eq. ( |l0| ) as well as, in an approximate, way, the core condi- 
tion via the truncated eq. (|ll| ) or ( |l3|) together with the approximation (^), we thus arrive at a set of equations 
implementing hrt with the loga/ ORPA-like closure (||) on the two-particle level well suited for numerical processing. 
While these expressions lend themselves to discretization in a straightforward way, it is computationally much more 
convenient to cast the pde in a form superficially resembling a quasi-linear one || so that an implicit finite-difference 
scheme requires only the inversion of a tridiagonal matrix. The re-writing we adopted — detailed in appendix |a|, 
very similar to the one of ^ — results in the introduction of an auxiliary function f(Q,g) so that the PDE implied 
by eqs. (|^) and (0) can be written in the form 

Q) = d oo[/, Q, g] + doi [/, Q, q] §- g f(Q, 0) + d 02 [f, Q, g] ^f(Q, q) ■ (14) 

Now it is easy to demonstrate the PDE's stiffness for subcritical temperatures: recalling the definitions of appendix [a|, 
we find that the inverse isothermal compressibility 1 /kt of the system with potential v can be written as 1/kt = —g 2 x 
u>(0, g)/e(0, g), where e(Q, g) + 1 = e(Q, g) oc exp(/(Q, g) Uq(Q, g)). For subcritical temperatures there is a density 
interval [g v , g{\ where 1/kt = Q (implying diverging E and hence /); here, g v and gi are the densities of the coexisting 
vapor and liquid, respectively (recall that hrt yields rigorously flat isotherms within the coexistence region, with 
binodal and spinodal coinciding in three dimensions [Q). By construction, however, the limit Q — > is a continuous 
one (cf. eqs. (||) and (0)) so that /, e and related quantities must be large already well before Q = is reached; 
at the same time, the RG mechanism introduced by the definition of via eq. (0) precludes any divergence at 
non-vanishing Q. Considering the region of the (Q, g)-plane where / and e are large, we easily find from the explicit 
expressions of appendix ^ that the d^i, and hence df(Q,g)/dQ, are of order e 1 ; restricting ourselves to a specific 
density g and sufficiently small Q (so that uo(Q, g) = 1; extension of the argument to a larger Q-range is cumbersome 
but straightforward) we can write 
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- e/W..)do(0) , 

dQ 

where do is now of order unity. Inspection of the solution of this ODE immediately shows that the average of do (Q) 
over the interval [Qi, Q 2 ], < Q\ < Q2, is rigorously bounded from above by exp(— /(Q 2 , g))/(Q2 — Qi) or else there 
were a singularity in that Q-interval; translating back to / we see that, while \df(Q, g)/dQ\ is still of order e, / must 
be a rapidly oscillating function of Q (with a period of order e _1 ), the average slope of which is much smaller, viz. of 
order 1/Q. 

It should be noted that this stiffness is not an artifact of the re- formulation of the PDE as summarized in appendix |A| 
but is manifest just the same when directly solving the PDE for the modified free energy A^\g) rather than that 
for the auxiliary f(Q,g) p^ ]. The above argument relies only on some general properties of hrt in the current 
formulation as applied to one-component fluids: the divergence of the isothermal compressibility in the coexistence 
region (the reproduction of which is one of HRT's main achievements), continuity of the limit Q — > 0, and the 
suppression of divergences for Q > as a result of the RG mechanism implemented via the truncated potential v^K 
The essential additional ingredient, viz. the behavior of the ratio of the Q- and the p-derivatives as the divergence 
in the compressibility builds up, while obvious in the formulation via /, is not easily seen in terms of A^\g); this, 
however, comes as no surprise since / is essentially the free energy's derivative with respect to Q so that we would 
have to reason about third- and second-order derivatives rather than first- and second-order ones if we were to repeat 
the arguments without resorting to the re-writing of appendix |X| 

III. OVERVIEW OF THE PROGRAM 

With the theory and approximations outlined in the previous section, we are now in a position to undertake the 
task of implementing hrt for one-component fluids in fully standards conforming Fortran-90; the only non-standard 
feature we make use of is the availability of the special values NaN and ±Inf for numerically undefined values and 
signed overflows, respectively, as defined in the floating-point standard IEEE 754. These requirements should not 
pose a serious restriction for our program's possible users: after all, Fortran-90 compilers have been available for a 
wide range of platforms for several years, and the desired floating-point behavior can usually be requested — albeit 
at a small performance penalty — via compiler switches. While our implementation is more appropriately described 
as a collection of mutually compatible building blocks rather than as a monolithic program so that the details of the 
numerical procedure are best left to these parts, for the combination of different selections to work all versions of all 
the modular constituents must adhere to a common model of the computation: 

Most obviously, we have to make the transition from the PDE's domain, viz. the infinite strip [0, 00) x [g m i n , ftnax], 
to a discrete mesh defined by a finite number of points in the (Q, p)-plane. Evidently, the placement of these "nodes" , 
as we shall call them, is of utmost importance for the quality of the discretization so that it is only natural to define 
e#, the central parameter governing all of the numerics, in terms of the properties of this mesh: the coarser a mesh 
we chose, the larger e# will be. — In principle, the locations of the nodes, the data structures of which are organized 
in linked lists, can be chosen freely; in particular, the cut-offs Q of all the systems in such a list of nodes are not taken 
to necessarily coincide, even though this is usually the case except for a low-density boundary at p — 0. As for the 
densities of the nodes, implementation of the core condition via the truncated eq. ([11]) and eq. (|12|) makes anything 
but constant (though not necessarily equispaced) density values impractical; if the grid is to be refined for low Q, 
additional nodes must be inserted at the same densities in all the node lists in the calculation. — After initialization of 
the nodes' data structures, solution of the PDE proceeds by applying a (possibly iterated) predictor-corrector scheme 
to generate an approximate solution for the nodes most advanced towards Q = from the information available 
through the node lists at higher Q; in the interest of the code's simplicity, the number of such node lists has been 
fixed to exactly three, thus facilitating determination of appropriate step sizes AQ. 

As mentioned before, the goals set for our implementation of hrt necessitate a fully modular architecture of our 
program; and while we did not want to forgo the well-known advantages of Fortran-90 for the numerical work, 
experience with prior versions of our code taught us that the kind of flexibility we need cannot be accommodated 
within the rather rigid framework Fortran-90's modules with their one-way flow of information provide. Instead we 
opted for a simple meta-language^ for self-configuring construction of code customized to the chosen combination of 
approximations and the physical system at hand, at the same time enhancing readability and maintainability of the 



2 Available on the world wide web from ittp: //purl . oclc . org/NET/arf g/ 
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source and encouraging modularization; for a more detailed discussion of this approach and the numerous technical 
advantages it affords we refer the reader to [f22f . 



A. Main parts 



Our software can best be understood as a collection of mutually compatible and freely exchangeable building blocks 
corresponding to the underlying physical and mathematical notions; the resulting natural organization of the code 
cleanly separating conceptually unrelated approximations is a direct consequence of our adoption of a meta-language 
and the use of automatic code generation techniques. The implementation's modular constituents, henceforth dubbed 
"main parts", must, however, be clearly distinguished from Fortran-90's modules: in general, there is no simple 
mapping from main parts to modules, and every main part may give rise to any number of modules, incorporating 
all the information available within the code base. 

In the following sub-sections we take a closer look at some of the main parts, their physical meaning, the algorithms 
and approximations implemented, and at some of the information they make available to the other parts; we will, 
however, exclude from this discussion the program's infrastructure, e. g. the implementation of logging, of reading 
and parsing of options files, handling of node lists and the definition of a versatile, lossless and storage-efficient 
albeit platform dependent file format for the results at Q — Qq. In a similar vein, we only mention the assortment 
of accompanying tools for reading these files and dumping their content in human-readable or Mathematica-usable 
form, for locating the critical point or calculating phase diagrams. — Thus only main parts "potential" , "reference" , 
"ansatz" and "solver" remain to be discussed: 



B. Properties of the potential 

First and foremost, we obviously have to provide information on the fluid's potential v — v lcf + w and its properties: 
this is the purpose of the main part labeled "potential" . Just as the full potential is a sum of a reference part w rcf 
and a perturbational part w, the functions and parameters to be provided by this main part fall into two distinct 
categories, pertaining to either v or w; in addition, as the temperature enters the calculation only as a pre-factor 
to w, viz. via (f> = —(3w, the inverse temperature (3 is also defined here. 

As far as the reference system is concerned, restriction to hard spheres (for the rationale cf. section ||) means that 
only a function returning the hard sphere diameter <j{q) and a flag indicating any deviation of o~(g) from the unit of 
length need to be made available. — A similar parameter pertaining to the perturbational part w of the potential, viz. 
a flag indicating any density-dependence of w, also plays an important role in many parts of the program as substantial 
simplifications and, in many cases, significant speed-ups by caching previous results are possible whenever 4>(Q, g) 
only depends on Q. In addition, at every cut-off Q the program must have access to the Fourier transforms w(Q, g) 

and 4>(Q, g) as well as the derivatives d(f>{Q, g)/dQ and d n (i(Q, g)j /dg n , whereas powers of the volume integral, 

<M0, £?)", and their derivatives d n (0(0, g)j /dg n obviously do not depend on Q (here, m and n are appropriate 

integers known during code construction). 

For the benefit of the PDE-solving algorithm, this main part also has to set a parameter A[„] related to the maximum 
relative curvature of the second Q-derivative of 4>{Q, g), defined in such a way as to coincide with A for the square 
well potential W SW [- £ . A ^I. Indeed, most of our efforts so far (l2) have concentrated on this particularly simple type of 
potential given by 



w sw[ - e ' x < a] (r) 



—e: r < \ a 
: r > Act , 



where the hard sphere diameter ct of the reference part and the strength e of the potential are usually chosen as units 
of length and energy, respectively. — Another type of potential we have implemented is a generalized step potential, 
i. e. a succession of stretches of constant w(r); more specifically, the perturbational part w(r) of an n-step potential 
of this type is just a sum of square-well potentials, w — Y^i=i w sw ^-~ ei ' Xi '' 7 \ \i < Xi+i, 1 < i < n; again, we have 
only considered the ^-independent case, while all the potential's parameters should be assumed functions of g when 
modeling a specific physical system. — We have also implemented a ^-independent hard-core Yukawa potential v hcy , 

^hcyf-eo ,-e,z,<r]^ = / - £ / : r < a 



- e <L e -z{r-*) ■ r> G ^ 
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where the parameter eo, the value of — w inside the core, defaults to e = — w hcy (a+), which again is usually chosen 
to coincide with the unit of energy; any mismatch between eo and e dominates w hcy (k) for large k and is found to 
render unstable at least the numerics. 

Regarding the stability of the pde, recall from section || that an attractive potential (so that w(0, g) < 0) is a 
necessary though not sufficient condition for the stability of the pde. 



C. Hard sphere reference system 

Due to the specialization of u rcf to hard spheres, the reference system enters the expressions of section |l| only through 
the direct correlation function c 1 ^ 1 , implementation of which is the task set for main part "reference": drawing upon 
the information from main part "potential", only some initialization code and functions for the evaluation of c| (Q, g) 
and dc T 2 ci (Q, g)/dQ have to be exported. 

In our program we have so far included two different versions implementing the Percus-Yevick p3[ approximation and 
the Henderson- Grundke |^4| description; note that all results reported here have been obtained using the Henderson- 
Grundke '■ in a theory relying on internal consistency conditions like eq. (Q) as heavily as hrt, the thermodynamic 
inconsistency present in the Percus-Yevick solution seems particularly undesirable. 



D. Discretization, boundary conditions, and other approximations 



Main part "ansatz" , where all the approximations on the physical and mathematical level are combined to jointly 
define a reasonable numerical model of hrt, is at the very core of the PDE-solving machinery: for the potential the 
perturbational and reference parts of which have been described in the two previous sub-sections, the hrt-pde is 
discretized and solved according to a given set of approximations and on the mesh defined by the node-lists served by 
main part "solver" (sub-section III E| ). More precisely, "ansatz" provides a set of facilities in the form of subroutines 
with standardized interfaces implementing the various stages of the computation, viz. initialization of the node lists 
at Q = Qoo and solution of the PDE according to a predictor-corrector full approximation scheme. Note, however, 
that the code must accommodate the possibilities of both iterating the corrector step (which may allow reaching the 
numerical quality indicated by e# with somewhat larger step sizes, thus speeding up the calculation) and of discarding 
part of the solution should e#-based criteria not be met; to aid "solver" in these decisions, care has to be taken to 
detect and signal numerical anomalies. Once a step's results have been accepted, "ansatz" may perform additional 
manipulations of the data structures; most importantly, the re-scaling of all quantities affected by exponentiation of 
/ necessary whenever / is large (cf. our discussion of the pde's stiffness in section |||) is adjusted only when the last 
corrector's result has been accepted. 

Due to the eminent role of the consistency condition eq. (Q) in constructing a closure to the underlying ODE (||), the 
pde ([l4|) for f(Q, g) is of first order in Q and of second order in g; assuming the lowest possible number of nodes in 
the discretization (extension to higher order is straightforward), we need at least a 2 x 3 set of nodes. According to the 
general model of the computation presented in section [II, however, we instead keep a third node list in order to allow 
monitoring of second Q-derivatives, so that we use a discretization on the 3x3 grid schematically presented in fig. [I], 
including information available via that additional node list. Locally, the discretization is derived from an expansion 
about the midpoint of the nodes labeled (22) and (32) in the schematic [l], evaluating the second ^-derivative along 
the line of constant Q through this point (thin horizontal line in fig. [l]) by estimating the data at the intersection 
with the lines of constant density by interpolants defined from node triples (il) and (i3), respectively; the resulting 
finite difference approximation is applied to every set of three adjacent node-triples, substituting suitable boundary 
conditions at g min and f? m ax- 

As indicated in fig. |], Q is not necessarily constant along a given node list, whereas the stability of the numerical 
scheme may impose certain geometrical constraints regarding the possible locations of the nodes, e. g. for ensuring that 
the Courant-Friedrichs-Lewy criterion |25}] is met or for maintaining convexity of the remaining integration region; 
a suitable representation of these constraints is exported and must be taken into account by main part "solver" . 
If the latter decides to insert nodes at intermediate densities, the code for initializing the inserted data structures 
and for interpolating appropriate quantities is negotiated between the main parts, depending upon the order of the 
interpolation formulas available. A further consequence of having non-constant Q is that some parts of the density 
range may reach Q w Qq earlier than others; in this case, the corresponding nodes are locked, preventing further 
modification, and all of the converged nodes except those necessary for providing a boundary condition for the 
remaining density interval are removed from the node lists available to main part "ansatz" . 
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In addition to the discretization of the hrt-PDE (|lj) discussed so far, the implementation of the core condition 
along the lines of section || and appendix [a| is also of interest. Relegating discussion of the choice of appropriate basis 
function u n , 1 < n < N cc , to appendix [b], we only point out the extremely slow convergence of the 2"<W-mtegrals @ 
that have to be evaluated at Q — Qoo', furthermore, as the integrand is temperature dependent for k > Q 00l these 
integrals have to be evaluated for every isotherm — a problem that might be sidestepped by adopting the original 
implementation's strategy of consistently using the results for Q — > oo rather than those valid at Qoo for initialization 
even though such an approach introduces a discontinuity at Q = Qoo- Also, with the usual choice of Qoo ~ 10 2 /er, 
integration merely up to k = Qoo can hardly be deemed sufficient; an appropriate upper integration limit can instead 
be found by comparing the integrand's asymptotic behavior with e#. — Main part "ansatz" also has to identify 
quantities suitable both for monitoring convergence of the full approximation scheme and for choosing appropriate 
step sizes AQ and Ag, and to make available code fragments for the inspection of nodes in various stages of the 
computation as well as a description of the boundary conditions at g m i n and g max including mandatory settings for 
either of these parameters if necessary; in particular, most implementations require g m i n = in order to be able to 
use the divergence of the ideal-gas term in as a Q-independent boundary condition for / (cf. appendix [A|) . 

It is this main part that defines the formulation of HRT and the set of approximations used in the calculations; of 
the numerous versions of this main part we produced only a few are to be mentioned here: both the re-implementation 
of the original program's approximations for the core and boundary conditions and the approach combining the PDE 
with = at all densities including g ma , x , while mathematically inconsistent, retain thermodynamic consistency 
at least in some approximate way (cf. our discussion of the decoupling assumption in section Q); in addition to these, 
the two possible approaches at least mathematically meaningful, viz. the odes directly following from decoupling and 
the pde resigning on the core condition for the benefit of the compressibility sum rule (Q) (with the loga/orpa 

prescription 7^ (g max ) — as high density boundary condition) will also be used in section |y|'s presentation of the 
results our software yields. 



E. Criteria for positioning of nodes 



If main part "ansatz" is to provide a discretization on whatever mesh is handed to it, it is the task set for "solver", 
the last of the main parts to be discussed here, to define this very mesh and to keep track of the numerical solution's 
quality. Based primarily upon the value of e# and respecting any restrictions exported by "ansatz", step sizes AQ 
and Ag have to be chosen and checked for compatibility with the solution generated, iterating or discarding steps if 
certain criteria are not met; whenever "ansatz" signals an exception the last step is discarded, accepting the data in 
the node list corresponding to labels (2i) in fig. [I] as the best approximation to the solution for Q — > 0. At the same 
time, care has to be taken to locate and identify any problems in the solution, i. e. parts of the (Q, p)-plane where the 
solution found does not appear smooth on the scales set by the step sizes, the most basic assumption underlying any 
finite-difference calculation; whenever this assumption no longer holds, the algorithm will react by locally reducing 



AQ and Ag, inserting node triples (cf. sub-section HID) in order to achieve the latter. — Once we find any nodes 



already holding the final results for their respective densities they must be taken care of as discussed in sub-section 



HID; integration of the PDE is ended when there is a node with Q < Qo for every density in the calculation, or when 



"ansatz" requests an end either because an error condition has occurred (v. s.) or because the current node list is 



sufficiently close to Q = Qo already. — As noted in section III, the intimate link between this main part's task and 
the numerical quality of the solution generated makes it natural to here define e#, the central parameter governing 
the numerics, and it is this part of the program that relies upon e# and the associated criteria the most; other main 
parts use e# for little more than for switching between full analytic expressions and asymptotic expansions (for a 
slightly atypical example of which cf . appendix |b|) . 

Of the two implementations of this main part, one has been written in the hopes of being able to avoid the 
problematic region of large /(Q, g) altogether, as is indeed possible for some similar PDEs. This implementation 
makes full use of e#, relying on numerous criteria to control the calculation; in the following discussion the notation 

Py^ refers to customization parameters that should usually be taken as real numbers of order unity. A case in point 
is the choice of the density grid: even though this is not necessary, we decided to always start with an equispaced set 
of g- values ranging from g m i n to g max _; the number N e + 1 of such density values is related to e# by 

* j (^?max £?min) 

8 = [q] ' 

reflecting the importance of second p-derivatives in the solution of the PDE as well as the static nature of this set 
of densities due to the I-approach to the core-condition. — Once e# has been determined, the system is ready to 







start determining appropriate step sizes AQ; in particular the assumption that the potential v(r, g) introduces length 
scales only in the ra nge fr om a{g) to X[ v ](g) cr(g), where A[„] is related to the fourth derivative of v(k, g) with respect 



to k (cf. sub-section [II B ), places an upper bound on the admissible step sizes, viz. 

AQ < 



1 o l A Q] 



X[ v ] a(g) 

On the other hand, for a finite difference scheme to be meaningful at least a certain number of bits must remain 
significant in evaluating the differences, which implies a lower bound on AQ proportional to Q, and the solution has 
to be smooth on the scales defined by the mesh, which also rules out abrupt changes in the step sizes; consequently, 
the ratio of two consecutive AQ steps at the same density is restricted to lie between p^^ and 1 /Pr»tS ■ ^ n a similar 
vein, considering smoothness in the ^-direction we have to postulate that (Q(22) + Q(32))/2 is greater than either of 
Q(3l) and Q(33), where the labels coincide with those of the nodes of fig. ||; this condition, unlike the other rules 
mentioned so far, does not limit the step sizes AQ at any density g but rather determines whether Ag should be 
reduced by the insertion of nodes at an additional density. But the most important criteria for choosing AQ come 
from monitoring the solution generated: for every monitored quantity x we make sure that 



dQ 7 



AQ<Je # pi AQ] 



IMIq 

IWIq = max fc>Q \x(k, g)\ ; 

the quantities taken for x are, of course, defined by "ansatz" , and a usual choice is x(Q, g) € |/(Q, g), 1/K^(Q, g)\. 

A different set of quantities y, also chosen by "ansatz" and usually comprising just y(Q,g) — f(Q,g), is used to 
monitor the convergence of the predictor-corrector scheme and to determine whether or not the corrector should 
be iterated: denoting the absolute difference of consecutive approximations of y divided by \\y\\Q by Ay, itera- 
tions are performed until Ay < e#p^y° nv \ and the ratio of two consecutive step sizes is bounded from above by 

(e # p{ / AQ1 ) 1/ ( 2+P ^ 1 ) /y/AWy, where A« y is Ay evaluated after the first corrector step. According to simple heuris- 
tic arguments regarding the convergence of corrector iterations and ignoring the effect of other criteria, an average of 
Pn-^ calls of the corrector can be expected to solve the difference equations to within e# , and a setting of P 1 ^^ > 1 
may significantly speed up the calculation by allowing larger steps to be taken without loss of accuracy. — After 
finding and tentatively using a candidate AQ, we still have to check that the assumptions leading to that particular 
choice for AQ actually hold; to this end we re-evaluate all the criteria with the obvious exception of the one involving 
A^y after the predictor and discard the step unless a slightly smaller step size, viz. AQ p]^®l ld (Pdf S ca r d ^ P asses 
the tests. If no step size can be found satisfying all the constraints, the calculation is terminated. 

While the above set of prescriptions for finding suitable node locations has proved indispensiblc in understanding 
the behavior of the PDE's solution, the oscillatory nature of /(Q, g) invariably linked to the build-up of the isothermal 
compressibility's divergence for sub-critical temperatures prevent its use for /3 > f3 c : considering even the modest 
value / ~ 10 3 , a AQ would have to be smaller than e~ 10 ~ 10 -430 , which is obviously completely useless for any 
practical implementation. Thus, even though it means loosing control over the level of accuracy in the solution, 
we have also implemented a version of this main part with predetermined step sizes that just happen to often be 
sufficient for reaching Q = Qo even well below the critical temperature while reproducing the overflow necessary for 
kt's divergence in a density interval the edges of which may then be identified with the coexisting phases' densities 
g v and gi . Recalling the behavior of / wherever it is large we obviously have to drastically reduce AQ as we approach 
Qo; for this we use the very prescription introduced by the authors of || and evidently underlying all later published 
hrt calculations. 

One last aspect of this main part to be mentioned regards the choice of Qoo in both implementations: As the only 
reasonable initial condition for the core condition assumes that the structure at Qoo is basically the same as that for 
Q = oo (so that c^°°^ — c£°^ — c' 2 cf or, equivalently, 7^°°'' = 7n°°' > = 0, n > 0) and the same set of parameters must 
also be used for the initialization of nodes at Q = Qoo + |AQ| (the nodes labled (li) in fig. [l] in the first step), it 
is preferable to have (g) / dQ — at Q = Qoo! from eq. ( pi] ) one immediately concludes that this is equivalent 

to w(Qoo) = whenever using the decoupling assumption. It is left to main part "ansatz" to decide whether Qoo 
should be determined in this way, thereby necessarily introducing ^-dependent Qoo when dealing with a p-dependent 
potential. 
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IV. DISCUSSION AND CONCLUSION 



In the preceding sections we had to introduce a number of approximations, some of which may seem rather less 
justified; their respective importance for and bearing on our program's predictions of structure and thermodynamics 
of simple liquids now remains to be assessed. As an exact solution with which to compare numerical results is 
lacking, for hrt as implemented by our software package to be considered a reliable tool well applicable to realistic 
physical potentials along the lines outlined so far it is necessary to demonstrate the limited effects variations in the 
numerical recipe have and to compare the results obtained with those available by other means for certain potentials. 
For this contribution we turn to the hard-core Yukawa fluid with z = 1.8/ct for illustration, a system that has also 
been considered in a recent study B comparing various thermodynamically consistent approaches including hrt to 
short-ranged potentials. 

Among the approximations just mentioned, an important one is the necessary truncation of eq. (Kw to a finite 



number N cc + 1 of basis functions u n and expansion coefficients 7n , n — 0, . . . ,N CC ; in fact, not only must N cc 
be finite, it should also be rather small if evaluation of the slowly-convergent I^-integrals at Q = Qoo is not to 
dominate program execution time. An obvious test for the minimum number of basis functions to keep is to look at 
the iV cc -dependence of the phase behavior predicted as summarized in table ||, where the inverse critical temperature 
fic, the critical density g c , and the coexisting densities at (3 — 0.9/e are recorded; the latter temperature is sufficiently 
far away from the critical point so that the differences in g v and gi are not merely to be attributed to the differences in 
/3 C while maintaining sufficient separation of the binodal from the boundaries at p m i n = and £ max {v. i). From this 
we find inclusion of the core condition to be of vital importance in determining the fluid's phase behavior, while non- 
negligible variation of the results remains even for the highest iV cc -values considered; however, the amount of variation 
especially in (3 C drops significantly as soon as iV cc exceeds 5. — The real test for the use of the truncated eq. ( p"l| ) 
and the ad hoc approximation ( |l2| ) (easily demonstrated to be inadequate at least for certain potentials Q), is the 
pair distribution function g(Q°)(r) as obtained from the final values of the expansion coefficients 7^ ' by performing 
the inverse Fourier transformation implied by the Ornstein-Zernike equation (|5|) : Employing the odes following from 
the consistent application of the decoupling assumption in order to isolate the approximations pertaining to the 
implementation of the core condition from other effects (v. i.), the g^°'(r) so obtained generally takes on rather large 
values for small r while remaining within a few percent of the contact value g^°\a-\-) for larger r up to a—; while 
increasing N cc usually does not considerably reduce the magnitude of g^°'(r) for r close to 0, it instead extends the 
r-range of rather small g(Q°\r) to ever smaller r. At high density there is no substantial improvement in g^°\r, g), 
r < a(g), for N cc > 5, nor for N cc > 7 at low density, a finding corroborated by direct inspection of the final values 

of the expansion coefficients "fi^ ^ |22[ . Accordingly, 7V CC should probably be chosen no less than 7 (corresponding to 
a 6 th order polynomial in r for C^-*(r) inside the core), whereas N cc = 5 may still be sufficiently accurate for some 
applications; for iV cc < 5, on the other hand, we cannot expect significantly better results than those from consistently 
solving the PDE without implementing the core condition (cf. sub-section HID), which, of course, runs much faster 
and at least does not rely on inconsistent approximations. But note that the core condition is always but poorly met, 
irrespective of N CCl whenever f(Qo, g) is large (corresponding to the coexistence region or the critical point's vicinity in 
implementations relying on a PDe). On the other hand, when solving the pde without implementing the core condition 
at all, g(®°'(r), r < a, can, of course, become arbitrarily large; e. g. for = 0.7/e and g = 0.9/a 3 , g(Q°>(r) = —3.26 
inside the core while the contact value is g^°\a+) = +1.91. All in all, while systematic shortcomings in the pair 
distribution function g^°\r) itself cannot be avoided in an implementation relying on eqs. jll|) and ( |l2| ) with finite 
N cc , we needs must keep the core condition in the calculation due to its bearing on the phase behavior predicted; this 
is somewhat at variance with earlier findings || indicating only a modest influence of the core condition upon the 
results, a finding expressly referred to in p6[ |. 

On the other hand, as pointed out already in section [n|, retaining the core condition is possible only when adopting 
the decoupling assumption, viz. cS®\g) = 0, an assumption that itself suffices to decouple the PDE to a set of odes 
(which is why we can use a^'(gmax) = as a boundary condition). But if such a procedure is to be considered 
harmless, the results of the consistent application of this approximation to the closure (|^) should not differ much from 
those of the pde applying the decoupling assumption only in order to get rid of the second X-integral in eq. ([ll]). 
However, the calculations summarized in fig. || clearly show that the two approaches yield very different results so 
that we cannot rule out a non-negligible effect on the structural and thermodynamic properties predicted: Most 
importantly, the odes cannot reproduce well-defined phase boundaries, and they even yield slightly negative inverse 
compressibility 1/kt in what would otherwise be the coexistence region. Preserving the structure of the pde, so that 
thermodynamic consistency is at least partly implemented by the pde's coefficients doi a s defined in eq. (Al), seems 
sufficient to remedy these deficiencies; at any rate, we have to accept the decoupling assumption as indispensable for 
the implementation of the core condition. 

Let us now turn to a brief discussion of several other aspects of hrt in its present formulation and the potential 
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problems they may present; as most of these effects are much more prominent for the square- well and multistep 



potentials defined in sub-section [II B , we relegate more thorough treatment to [|12|,|22j . — Clearly, the importance of 
retaining the structure of a PDE mandates closer examination of the properties of the set of densities, especially since 
the terms corresponding to the finite difference approximation to the operator (d 2 /dg 2 ) have been found the primary 



limiting factor for the quality of the numerical solution (cf. the definition of e# in sub-section HIE); while obvious for 
subcritical temperatures at low Q where the near-discontinuity of / at the phase boundary betrays the smoothness 
assumptions underlying finite difference schemes, this is true even for rather small f(Q, g). For our hard-core Yukawa 
system, however, the results' stability with respect to a variation of the density grid or the location and nature of the 
high density boundary (cf. table |n|) is rather satisfactory as long as Q max — Qi exceeds several Ag: just as expected 
from fig. ^, the ode used at the boundary forces the corresponding density to lie outside the coexistence region, which 
carries over to nearby densities by virtue of the pde's discretization. Despite the identical phase behavior found for 
different boundary conditions as evidenced by table |j, the question of which condition to impose at g ma , x is far from 
irrelevant; indeed, when imposing the LOGA/ORPA-condition 7q = at p max without making use of decoupling at 
all, inspection of the solution generated close to g max clearly shows the need to replace the original implementation's 
approach by an ansatz less inconsistent. — Similarly, the stiffness of the pde for subcritical temperatures shortly 
touched upon in section || is easily detected by employing the criteria discussed in sub-section III E| ; still, for u hcy 



with z = 1.8/(7 and N cc = 7, using pre-determined step sizes and resigning on any control of the predictor-corrector 
scheme's convergence we obtain an approximate solution at Qq sufficiently stable outside the coexistence region; thus, 
while the solution generated in the region of large f(Q,g) necessarily differs from the true solution in a fundamental 
way, the influence of which on the data produced outside the coexistence region is not assessed easily, stiffness here 
appears less pressing a concern than for other systems Jl2| . 

In addition to the internal consistency of the results, we also have to make contact with data available by other 
means; as this paper's purpose is presentation of our software and discussion of some general aspects of hrt rather 
than a comprehensive study of the hard-core Yukawa system, we here restrict ourselves to the data of table I of 
H containing the predictions of various thermodynamically self-consistent liquid state theories including HRT in the 
original implementation, which are found to yield T c ranging from 1.193e/fcs to 1.219e/fcs (/3 C = l/fcsT c between 
0.820/e and 0.838/e), as well as the Monte Carlo (mc) result T c = 1.212(2)e/fc B (/3 C = 0.825(2)/e). As is apparent 
from our table |, our implementation's predictions for 1 < N cc < 4 fall precisely into this range and, for N cc £ {3,4}, 
are well compatible with simulation results; that very same iV cc -range, on the other hand, we have seen (v. s.) is 
characterized by gross violation of the core condition due to an insufficient number of basis functions retained in 
the truncated eq. (p"T|). As we further increase A cc so that the core condition is obeyed to a certain extent (v. s.), 
however, (} c drops dramatically to values far outside the range quoted in Jj|; while the trend of increasing (3 C evident 
for N cc > 7 indicates that hrt might match the MC predictions for A cc ~ 15, we have not performed these CPU 
intensive calculations. 

With this we conclude our superficial sketch of the software we have written and the appraisal of the results it 
typically produces as illustrated for the hard-core Yukawa potential v hcy with z = 1.8/er: it should be apparent that 
HRT in its current formulation, while presenting substantial difficulties discussed here as well as in [[12|22), is well 
capable of predicting structural and thermodynamic properties of simple one-component fluids; at the same time, 
the computational difficulties mentioned and the approximations introduced to render the numerics tractable cannot 
always be shown to be harmless so that great care has to be exercised in the interpretation of isolated results. The 
fully modular design of our program and the high degree of flexibility brought about by the adoption of a meta 
language and code construction techniques are key factors in facilitating additional evaluations, providing a means of 
separating different approximations' effects; at the same time, our implementation makes for a versatile tool for the 
systematic exploration of hrt for one-component fluids in its present or alternative formulations. 
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APPENDIX A: RE-FORMULATION IN NOT-QUITE QUASI-LINEAR FORM 

As noted in section |l[ it is advantageous to replace the highly non-linear pde in the modified free energy 
implied by eqs. (o), (pi), and the core condition, with a formulation akin to a quasi- linear one. In complete analogy 
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to [@], we define the auxiliary function f(Q, g) via 



in the ideal gas limit, i. e. for g — > 0, we immediately find from the divergence of the terms —1/g in C2° f and C^' 
that f(Q,0) — df(Q,0)/dg = 0, which is a convenient boundary condition most implementations rely on. In the 
above definition we have taken advantage of some freedom regarding the /-term to reduce the number of floating-point 
multiplications; as Uq is usually chosen to be strictly proportional to 0, our / differs only by a constant factor from the 
choice of || which can be recovered by replacing the factor Uq(Q, g) by 4> 2 (Q, g); note, however, that the expressions 
we give here remain valid for a slightly more general choice of uq, allowing the proportionality of the basis function 
and the perturbational part of the potential to hold only up to order 2 — a freedom that might be exploited to use 
more appropriate basis functions, giving 

C^) a i ar 

ger range in r-space. Also we should point out that the condition 
of non-singular coefficients in the pde ([lij), the very reason for the restriction on the relation between uq and 
just mentioned as well as for the introduction of the term involving (j)/K^\ only fixes a minimum exponent for the 
Mo-factor with which to multiply / . 

Inserting the above definition for / into the relevant equations of section |n| and eliminating the expansion coefficient 
7o via the consistency condition (Q), we can easily re-cast the original pde in the form of eq. (14). Dropping the 
obvious arguments and with the shorthand notations 

£=1-^ = e /«o+^ s = e-l 
= O = 0(0, e ) Go = G (Q) (0, g) 

the coefficients of eq. can be written as E3] 
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When evaluating these expressions care has to be taken whenever (j>(Q, g) is close to zero as all terms of orders 
\/4> 2 {Q,g) and l/(f>(Q,g) must cancel; noting that (£ 2 0g/e0 3 ) - (0o/(X^ Q) ) 2 0), the coefficient of dK.^/dQ in d o, 
can be written in terms of the (90/9(5)-coeflicient, in our numerical work the calculation of the terms affected 
by the cancellation for small 4>{Q,g) proceeds via application of a fifth-order Taylor expansion of {4>q/K,^ 2 ) — 
(K.^ e 2 4>q / e cj) 4 ) — (2//0); one should note that, even though the criterion for switching between the full analytic 
expressions and said expansion depends on e#, the order of this expansion (and a similar one for e/0) is not increased 
for very small values of e#, which is one of the few hard-coded limitations of our program. 

From the given expressions for the doi two more aspects are obvious: (i) there is substantial further simplification 
for a ^-independent potential (as uq depends upon the density only through the perturbational part of the potential, 
and the basis functions entering Q^> only through the reference system's hard-core radius cr(g)), and (ii) the pde 
does not fall into the class of quasi-linear pdes due to the presence of the second p-derivative of x<f, in the coefficient 



APPENDIX B: BASIS FUNCTIONS FOR THE CORE CONDITION 

While the basis function Uo(r, g) oc w(r, g) in the closure (^) is fixed by the potential used, in principle there is ample 
freedom in choosing the u n (r), n > 1. Of course, when truncating eq. ( p"T| ) after 1 + N cc basis functions it is natural to 
ask for the set {u± (r ) , . . . , m/v cc (f)} to span the space of polynomials of order up to iV cc — 1 , so that u n (r) will generally 
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be a polynomial or order n — 1 in r; but whereas different polynomials of this type do not alter the function space, 
their choice has implications for the numerical properties of the matrix equations implementing the core condi tion a s 
well as, to a certain extent, for the convergence of the i"-integrals to be evaluated at Q — Qoo (cf. sub-section [II D| ), 
Other than the original implementation that relied on an affine transformation of the Legendre polynomials, we found 
it convenient to choose u n (r, g) simply proportional to a power of r and normalized to u n (0, g) = 1. Thus, dropping 
the obvious argument g, we have 



u n {r) 



n+2 fr\ 
47r<T a \a) 



71—1 



the Fourier transform of which is 



l n (k) = 



(ak) 



n+2 



7J7T x - 



3=0 



(n - j) 



r < a 
r > a , 



(ak)""- 3 cos ( o-fc + y 



an expression used for ak > n only for numerical reasons. For smaller k we rely on the expansion 

(-iy 



In(fc) = (n + 2) ^ 



3=0 



(n + 2j + 2) (2j + l)\ 



(ak) 



truncating the series after N n terms, where N n is the smallest number such that 



,2Af„+2 



n + 2N n + 4 (2iV„+3)! 
with a customization factor p^™' of order unity. 



[1] A. Parola, L. Reatto, Adv. Phys. 44 (1995) 211. 
[2] A. Parola, L. Reatto, Phys. Rev. Lett. 53 (1984) 2417. 
[3] A. Parola, L. Reatto, Phys. Rev. A 31 (1985) 3309. 
[4] A. Parola, D. Pini, L. Reatto, Phys. Rev. E 48 (1993) 3321. 
[5] A. Parola, J. Phys. C: Solid State Phys. 19 (1986) 5071. 
[6] A. Meroni, A. Parola, L. Reatto, Phys. Rev. A 42 (1990) 6104. 
[7] A. Meroni, L. Reatto, M. Tau, Mol. Phys. 80 (1993) 977. 
[8] M. Tau, A. Parola, D. Pini, L. Reatto, Phys. Rev. E 52 (1995) 2644. 
[9] C. Caccamo, G. Pellicane, D. Costa, D. Pini, G. Stall, Phys. Rev. E 60 (1999) 5533. 
[10] L. Reatto, A. Parola, J. Phys.: Condens. Matter 8 (1996) 9221. 

[11] F. Barocchi, P. Chieux, R. Fontana, R. Magli, A. Meroni, A. Parola, L. Reatto, M. Tau, J. Phys.: Condens. Matter 9 

(1997) 8849. 
[12] A. Reiner, G. Kahl, to be published. 
[13] L. Reatto, Phys. Lett. 72 A (1979) 120. 
[14] D. Chandler, J. D. Weeks, Phys. Rev. Lett. 25 (1970) 149. 
[15] D. Chandler, J. D. Weeks, H. C. Andersen, J. Chem. Phys. 54 (1971) 5237. 
[16] H. C. Andersen, J. D. Weeks, D. Chandler, Phys. Rev. A 4 (1971) 1597. 
[17] A. Parola, A. Meroni, L. Reatto, Phys. Rev. Lett. 62 (1989) 2981. 
[18] L. S. Ornstein, F. Zernike, Proc. Akad. Sci. (Amsterdam) 17 (1914) 793. 
[19] G. Stell, Phys. Rev. 184 (1969) 135. 
[20] G. Stell, J. Chem. Phys. 55 (1971) 1485. 

[21] H. C. Andersen, D. Chandler, J. Chem. Phys. 55 (1971) 1497. 

[22] A. Reiner, PhD Thesis, Technische Univers itat Wien, 2002. Available on the world wide web from 

http : //purl . oclc . org/NET/a-reiner/dr-thesis/ 
[23] J. L. Lebowitz, Phys. Rev. 133 (1964) A895. 



14 



[24] D. Henderson, E. W. Grundke, J. Chem. Phys. 63 (1975) 601. 

[25] R. Courant, K. O. Friedrichs, H. Lewy, Mathematische Annalen 100 (1928) 32. 

[26] D. Pini, A. Parola, L. Reatto, J. Stat. Phys. 100 (2000) 13. 



iVcc 


Pc e 


g v (/3* = 0.9/e) a 3 


Qi{P* = 0.9/e) a 3 


Qc o 3 




n (\a( ^q^ 


nii ^ 

V. -L -LOyO ) 






1 


0.82070(39) 


0.105(5) 


0.575(5) 


0.315(30) 


2 


0.82148(39) 


0.105(5) 


0.565(5) 


0.315(10) 


3 


0.82227(39) 


0.105(5) 


0.565(5) 


0.315(30) 


4 


0.82227(39) 


0.105(5) 


0.565(5) 


0.315(30) 


5 


0.77676(20) 


0.075(5) 


0.645(5) 


0.320(15) 


6 


0.75527(20) 


0.055(5) 


0.685(5) 


0.325(30) 


7 


0.75957(20) 


0.065(5) 


0.675(5) 


0.330(25) 


8 


0.77305(39) 


0.065(5) 


0.645(5) 


0.320(35) 


9 


0.78555(39) 


0.075(5) 


0.615(5) 


0.315(20) 



TABLE I. Dependence of inverse critical temperature p c = 1/ks T c , coexisting densities g v and Qi at j3 = 0.9/e, and critical 
density g c for a hard-core Yukawa potential with z = 1.8/a on the number of basis functions. The results reported have been 
obtained from pdes retaining 7V CC + 1 basis functions or (first line) not implementing the core condition at all, with e# = 10 -2 
1.0/<t 3 . We have checked that the differences summarized here cannot be explained by the 7V cc -dependence of the 



and £ m 

upper integration limits in evaluating the X ( - < ^° c '-integrals (cf. sub-section HID) 



method 


A^cc 


bound, cond. 


Pet 


q v (P = 0.9/e) a 3 


gi{p = 0.9/e) a 3 


g c a 3 


no core condition 




7o WJ (f? max ) =0 


0.831573(91) 


0.115(5) 


0.565(5) 


0.330(20) 


decoupling 


7 


lo^^ftnax) =0 


0.759429(91) 


0.055(5) 


0.675(5) 


0.325(15) 


decoupling 


7 


a (Q) ( ema x) = 


0.759429(91) 


0.055(5) 


0.675(5) 


0.325(15) 



TABLE II. Inverse critical temperature p c = 1/fcs T c , coexisting densities g v and gi at j3 — 0.9/e, and critical density g c for 
a hard-core Yukawa potential with z — 1.8/a as predicted by various combinations of approximations and boundary conditions 
at £> max = 1/cr 3 . Again, the results reported have been obtained from pdes retaining iV cc + 1 basis functions or (first line) not 
implementing the core condition at all, with e# = 10 -2 and £> max = 1.0/cr 3 . 
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p 

FIG. 1. Schematic of the grid used in the discretization or the pde (cf. sub-section HID). Assuming use of the three-point 
approximation for the second derivatives in the ^-direction, the discretization is generated from an expansion around the 
intersection of the thin horizontal line with the line of constant density joining the nodes labeled (i2). — According to the 
general model of the computation discussed in section a node list's Q-values may be ^-dependent, whereas the ^-values 
must coincide in all three node lists, though they need not be equispaced. 




FIG. 2. Comparison of inverse compressibility of the hard-core Yukawa system with z — 1.8/ a as obtained from the odes 
following from the decoupling assumption a^\g) = (thin lines) and from a pde (N cc = 7) inconsistently applying this 
approximation to the evolution of the core condition expansion coefficients 7^' only (thick line). In both cases the equations 
have been solved on an equispaced density grid with Ag = 0.1/cr 3 , and the decoupling assumption also served as a boundary 
conditon for the pde at £> ma x = 1.0/ct 3 . The critical temperature obtained from the pde and used in the middle plot is 
(3 C = 0.759497(24)/e. 



1G 



